If we treat it as “stationary” series, we can try to fit ARMA(p,q) series.
## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone 'zone/tz/2018c.
## 1.0/zoneinfo/America/New_York'
source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/lake.csv")
Lake <- ts(D, start=1875, freq=1)
length(Lake)## [1] 97
##
## 'tseries' version: 0.10-42
##
## 'tseries' is a package for time series analysis and
## computational finance.
##
## See 'library(help="tseries")' for details.
## Warning in kpss.test(A): p-value smaller than printed p-value
## KPSS ADF PP
## p-val: 0.01 0.24 0.025
## Series: Lake
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.7665 0.3393 9.1290
## s.e. 0.0773 0.1123 0.3861
##
## sigma^2 estimated as 0.4784: log likelihood=-101.09
## AIC=210.18 AICc=210.62 BIC=220.48
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.963 0.952 0.934 0.567 0.641 0.89 0.684
Analysis 1 (direct fit)
auto.arima() chooses AR(2) with min AICC.
AR(2) with constant mean was fit directly to data. With observation \(Y_t\),
\[
Y_t \hspace3mm = \hspace3mm \mu + X_t \\\\
X_t \hspace3mm = \hspace3mm \phi_1 X_{t-1} + \phi_2 X_{t-2} + e_t
\]
#--- 10-step prediction ---
Y <- Lake # Y should be TS object
Fit <- Fit1 # Fit you are using for forecasting
Y.h <- predict(Fit, n.ahead=10) # Y.h$pred is TS object
Yhat <- ts(cbind(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se), start=start(Y.h$pred), freq=frequency(Y))
ts.plot(Y, Yhat, type="o", col=c("black","red","blue","blue"))
abline(h=0)
abline(h=mean(Y), col="blue") ##- 1. Direct ARMA Fit
#--- Rolling 1-step prediction ---
Y <- Lake
Window.size <- 60 # must be less than length(Y)
ARIMA.order <- c(0,0,2) # model to fit each step
Yhat <- matrix(0, (length(Y)-Window.size), 3) # initialize what needs to be saved
for (i in 1:(length(Y)-Window.size)) {
Fit <- Arima( Y[i:(i+Window.size-1)], order=ARIMA.order) # <--- Fit ARMA(p,q)
Y.h <- predict(Fit, n.ahead=1) # one step prediction
Yhat[i,] <- c(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se)
}
Yhat <- ts(Yhat, start=time(Y)[Window.size+1], freq=frequency(Y))
#- Calculate prediction performance
Pred.error = Y[(Window.size+1):length(Y)] - Yhat[,1]
Pred.rMSE = sqrt( mean( (Pred.error)^2 ) ) #- prediction root Mean Squared Error
Pred <- matrix(c(Pred.rMSE, 1.96*sqrt(Fit$sigma2), mean(Pred.error)), 1, 3)
colnames(Pred) <- c("Pred.rMSE", "95%PI", "Mean")
Pred
#- Plot the prediction result with original data ---
ts.plot(Y, Yhat, type="o", col=c("black", "blue", "red", "red"))
#- Plot of the prediction error with ACF of pred error
layout(matrix(1:2, 1, 2)) # two plots side by side
plot( Pred.error, type="o", ylim=c(-Pred[2], Pred[2])*1.2, main="Prediction Error (Blue-Black)" )
abline(h=0)
abline(h=c(-1.96, 1.96)*sqrt(Fit$sigma2), col="blue", lty=2)
acf(Pred.error)
layout(1) # reset the layout## Pred.rMSE 95%PI Mean
## [1,] 0.8029921 1.526601 -0.07798477
If we choose that the level is going down linearly, we can try to fit a line with ARMA errors.
##
## Call:
## lm(formula = Lake ~ time(Lake))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.50919 -0.74760 -0.01556 0.75966 2.53409
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.278141 7.922614 6.977 4.02e-10 ***
## time(Lake) -0.024071 0.004119 -5.843 7.16e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.136 on 95 degrees of freedom
## Multiple R-squared: 0.2644, Adjusted R-squared: 0.2566
## F-statistic: 34.14 on 1 and 95 DF, p-value: 7.165e-08
## Series: Reg2$residuals
## ARIMA(1,0,1) with zero mean
##
## Coefficients:
## ar1 ma1
## 0.6671 0.3827
## s.e. 0.0937 0.1135
##
## sigma^2 estimated as 0.452: log likelihood=-98.72
## AIC=203.44 AICc=203.7 BIC=211.16
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.985 0.974 0.969 0.19 0.189 0.782 0.669
## Series: Lake
## Regression with ARIMA(1,0,1) errors
##
## Coefficients:
## ar1 ma1 intercept xreg
## 0.6682 0.3817 55.5443 -0.0242
## s.e. 0.0936 0.1136 18.0324 0.0094
##
## sigma^2 estimated as 0.4612: log likelihood=-98.66
## AIC=207.33 AICc=207.98 BIC=220.2
Analysis 2 (linear trend)
Linear trend was fit to \(Y_t\), then the residuals were fitted with ARMA(1,1).
Y_t = a+b t + X_t X_t
In other words, \[ Y_t \hspace3mm = \hspace3mm a + bt + X_t \\\\ X_t \hspace3mm = \hspace3mm \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} \]
#- forecast ARIMA
plot(forecast(Fit2b, xreg=1972:1982))
abline(a=Fit2b$coef["intercept"], b=Fit2b$coef["xreg"], col="red") #--- 10-step prediction ---
Y <- Lake # Y should be TS object
Fit <- Fit2b # Fit you are using for forecasting
Y.h <- predict(Fit, n.ahead=10, newxreg=c(1:10)+max(time(Y)) ) # Y have to have freq=1
Yhat <- ts(cbind(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se), start=start(Y.h$pred), freq=frequency(Y))
ts.plot(Y, Yhat, type="o", col=c("black","red","blue","blue"))
abline(h=0)
abline(a=Fit$coef["intercept"], b=Fit$coef["xreg"], col="red") ##- 2. ARMA with Linear Fit
#--- Rolling 1-step prediction ---
Y <- Lake
Window.size <- 60 # must be less than length(Y)
ARIMA.order <- c(1,0,1) # model to fit each step
Yhat <- matrix(0, (length(Y)-Window.size), 3) # initialize what needs to be saved
for (i in 1:(length(Y)-Window.size)) {
#- Fit ARMA(p,q) WITH Reg Line
Fit <- Arima( Y[i:(i+Window.size-1)], order=ARIMA.order, xreg=time(Y[i:(i+Window.size-1)]))
#- 1-step Prediction for next Year
Y.h <- predict(Fit, n.ahead=1, newxreg=1+max(time(Y[i:(i+Window.size-1)])) )
Yhat[i,] <- c(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se)
}
Yhat <- ts(Yhat, start=time(Y)[Window.size+1], freq=frequency(Y))
#- Calculate prediction performance
Pred.error = Y[(Window.size+1):length(Y)] - Yhat[,1]
Pred.rMSE = sqrt( mean( (Pred.error)^2 ) ) #- prediction root Mean Squared Error
Pred <- matrix(c(Pred.rMSE, 1.96*sqrt(Fit$sigma2), mean(Pred.error)), 1, 3)
colnames(Pred) <- c("Pred.rMSE", "95%PI", "Mean")
Pred
#- Plot the prediction result with original data ---
ts.plot(Y, Yhat, type="o", col=c("black", "blue", "red", "red"))
#- Plot of the prediction error with ACF of pred error
layout(matrix(1:2, 1, 2)) # two plots side by side
plot( Pred.error, type="o", ylim=c(-Pred[2], Pred[2])*1.2, main="Prediction Error (Blue-Black)" )
abline(h=0)
abline(h=c(-1.96, 1.96)*sqrt(Fit$sigma2), col="blue", lty=2)
acf(Pred.error)
layout(1) # reset the layout## Pred.rMSE 95%PI Mean
## [1,] 0.7817429 1.497827 0.2225482
Let’s look at the difference between the observations: \[ diff(Y) = Y_t - Y_{t-1} \\ \\ \]
## Warning in adf.test(A): p-value smaller than printed p-value
## Warning in pp.test(A): p-value smaller than printed p-value
## Warning in kpss.test(A): p-value greater than printed p-value
## KPSS ADF PP
## p-val: 0.1 0.01 0.01
## Series: Lake
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## 0.6385 -0.5349 -0.3514
## s.e. 0.1345 0.1445 0.1055
##
## sigma^2 estimated as 0.4812: log likelihood=-99.88
## AIC=207.76 AICc=208.2 BIC=218.02
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.979 0.956 0.942 0.524 0.495 0.577 0.677
Analysis 3 (take difference)
Difference of \(Y_t\) was taken. The difference seems to be WN.
\[ \bigtriangledown Y_t \mbox{ is WN} \\ Y_t \mbox{ is ARIMA(0,1,0)} \]
#--- 10-step prediction ---
Y <- Lake # Y should be TS object
Fit <- Fit3 # Fit you are using for forecasting
Y.h <- predict(Fit, n.ahead=10) # Y.h$pred is TS object
Yhat <- ts(cbind(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se), start=start(Y.h$pred), freq=frequency(Y))
ts.plot(Y, Yhat, type="o", col=c("black","red","blue","blue"))
abline(h=0)
abline(h=mean(Y), col="blue") ##- 3. ARIMA Fit (Random Trend)
#--- Rolling 1-step prediction ---
Y <- Lake
Window.size <- 60 # must be less than length(Y)
ARIMA.order <- c(1,1,2) # model to fit each step
Yhat <- matrix(0, (length(Y)-Window.size), 3) # initialize what needs to be saved
for (i in 1:(length(Y)-Window.size)) {
Fit <- Arima( Y[i:(i+Window.size-1)], order=ARIMA.order) # <--- Fit ARMA(p,q)
Y.h <- predict(Fit, n.ahead=1) # one step prediction
Yhat[i,] <- c(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se)
}
Yhat <- ts(Yhat, start=time(Y)[Window.size+1], freq=frequency(Y))
#- Calculate prediction performance
Pred.error = Y[(Window.size+1):length(Y)] - Yhat[,1]
Pred.rMSE = sqrt( mean( (Pred.error)^2 ) ) #- prediction root Mean Squared Error
Pred <- matrix(c(Pred.rMSE, 1.96*sqrt(Fit$sigma2), mean(Pred.error)), 1, 3)
colnames(Pred) <- c("Pred.rMSE", "95%PI", "Mean")
Pred
#- Plot the prediction result with original data ---
ts.plot(Y, Yhat, type="o", col=c("black", "blue", "red", "red"))
#- Plot of the prediction error with ACF of pred error
layout(matrix(1:2, 1, 2)) # two plots side by side
plot( Pred.error, type="o", ylim=c(-Pred[2], Pred[2])*1.2, main="Prediction Error (Blue-Black)" )
abline(h=0)
abline(h=c(-1.96, 1.96)*sqrt(Fit$sigma2), col="blue", lty=2)
acf(Pred.error)
layout(1) # reset the layout## Pred.rMSE 95%PI Mean
## [1,] 0.7427918 1.498318 0.09050329
From Model1, 2, 3, estimate for \(\sigma^2\) was \(0.4812, \, 0.4612, \, 0.4784\) respectively.
Taking square root, that means \(\hat \sigma\) are \(0.6936, \, 0.6791, \, 0.6916\)
On the other hand, rolling 1-step prediction rMSE for the three models were \(0.8029, \, 0.7817, \, 0.7427\)
Why do they not agree? Which model is fitting the best?
What is over-fittting?